Milo: liver cirrhosis analysis

Introduction

In this notebook we demonstrate how to use Milo to detect abherrant cell states in diseased tissues, using a dataset of hepatic non-parenchymal cells isolated from 5 healthy and 5 cirrhotic human livers. Ramachandran et al. 2019 (GEO accessiion: GSE136103).

suppressPackageStartupMessages({
  library(tidyverse)
  library(irlba)
  library(DropletUtils)
  library(scater)
  library(scran)
  library(Seurat) ## just 4 loading the object
  library(miloR)
  library(SingleCellExperiment)
  library(patchwork)
  library(igraph)
  library(RColorBrewer)
  library(cowplot)
  })

Load data

We downloaded the dataset and annotations stored in Seurat object from here, as indicated by the authors.

load("/nfs/team205/ed6/data/Ramachandran2019_liver/tissue.rdata")

## Convert to SingleCellExperiment
liver_sce <- SingleCellExperiment(assay = list(counts=tissue@raw.data, logcounts=tissue@data),
                                  colData = tissue@meta.data)

liver_sce

class: SingleCellExperiment 
dim: 23498 58358 
metadata(0):
assays(2): counts logcounts
rownames(23498): FO538757.2 AP006222.2 ... CTA-126B4.7 LINC01423
rowData names(0):
colnames(58358): Healthy1_Cd45+_AAACCTGCAGTATCTG
  Healthy1_Cd45+_AACTGGTTCATGGTCA ...
  Cirrhotic3_Cd45-_TTTGTCATCCAGGGCT
  Cirrhotic3_Cd45-_TCTGGAAGTCATCCCT
colData names(10): nGene nUMI ... annotation_indepth
  annotation_lineage
reducedDimNames(0):
spikeNames(0):
altExpNames(0):

Preprocessing

Feature selection

dec_liver <- modelGeneVar(liver_sce)

fit_liver <- metadata(dec_liver)
plot(fit_liver$mean, fit_liver$var, xlab="Mean of log-expression",
    ylab="Variance of log-expression")

hvgs <- getTopHVGs(dec_liver, n=3000)

Dimensionality reduction with PCA

liver_sce <- runPCA(liver_sce, subset_row=hvgs, ncomponents=11)

plotPCA(liver_sce, colour_by="condition", ncomponents=3)

liver_sce <- runUMAP(liver_sce, dimred="PCA", ncomponents=2)

scater::plotUMAP(liver_sce, colour_by="condition", point_alpha=1,  point_size=0.5)

scater::plotUMAP(liver_sce, colour_by="dataset", point_alpha=0.3,  point_size=0.5)

scater::plotUMAP(liver_sce, colour_by="annotation_lineage", point_alpha=0.3,  point_size=0.5, text_by='annotation_lineage')

# scater::plotUMAP(liver_sce, colour_by='annotation_indepth', point_alpha=0.3,  point_size=0.5, text_by='annotation_indepth')

Notably, this dataset doesn’t appear to display a batch effect

DA analysis with Milo

We test for differential abundance between healthy and cirrhotic livers. We start by defining neighbourhoods with refined sampling on the KNN graph. We inspect the size of neighbourhoods.

liver_milo <- Milo(liver_sce)

## Build KNN graph
liver_milo <- buildGraph(liver_milo, d = 11, k=30)

Constructing kNN graph with k:30

## Compute neighbourhoods with refined sampling
liver_milo <- makeNhoods(liver_milo, k=30, d=11, prop = 0.05, refined=TRUE)

Checking valid object

plotNhoodSizeHist(liver_milo, bins=150)

Then we make a design matrix for the differential test, assigning samples to conditions.

liver_meta <- as.tibble(colData(liver_milo)[,c("dataset","condition")])

`as.tibble()` is deprecated as of tibble 2.0.0.
Please use `as_tibble()` instead.
The signature and semantics have changed, see `?as_tibble`.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

liver_meta <- distinct(liver_meta) %>%
  mutate(condition=factor(condition, levels=c("Uninjured", "Cirrhotic"))) %>%
  column_to_rownames("dataset")

Now we can count cells in neighbourhoods and perform the DA test.

liver_milo <- countCells(liver_milo, samples = "dataset", meta.data = data.frame(colData(liver_milo)[,c("dataset","condition")]) )

Checking meta.data validity
Setting up matrix with 2717 neighbourhoods
Counting cells in neighbourhoods

liver_milo <- calcNhoodDistance(liver_milo, d=11)

milo_res <- testNhoods(liver_milo, design = ~ condition, design.df = liver_meta, fdr.weighting = "k-distance")

Performing spatial FDR correction withk-distance weighting

Exploration of DA results

We can start by looking at some basic stats

pval_hist <- milo_res %>%
  ggplot(aes(PValue)) +
  geom_histogram(bins=50) +
  theme_bw(base_size=14)

volcano <-
  milo_res %>%
  ggplot(aes(logFC, -log10(SpatialFDR))) +
  geom_point(size=0.4, alpha=0.2) +
  geom_hline(yintercept = -log10(0.1)) +
  xlab("log-Fold Change") +
  theme_bw(base_size=14)

pval_hist + volcano

The distribution of P-values looks sensible and from the volcano plot we can see that we have identified some DA neighbourhoods at 10% FDR. We can visualize DA neighbourhoods building an abstracted graph

liver_milo <- buildNhoodGraph(liver_milo)

Calculating nhood adjacency

\

Error: unexpected input in "\"

## Save milo object and results
saveRDS(liver_milo,"~/liver_milo_20201008.RDS")
write_csv(milo_res,"/nfs/team205/ed6/data/Ramachandran2019_liver/liver_results_20201008.csv")

liver_milo <- readRDS("~/liver_milo_20201008.RDS")
milo_res <- read_csv("/nfs/team205/ed6/data/Ramachandran2019_liver/liver_results_20201008.csv")

## Load hvgs 
hvgs <- scan("~/liver_milo_hvgs.txt", "")

Making figures for the manuscript

colourCount = length(unique(liver_milo$annotation_lineage))
getPalette = colorRampPalette(brewer.pal(9, "Set2"))

umap_df <- data.frame(reducedDim(liver_milo, "UMAP"))
colnames(umap_df) <- c("UMAP_1", "UMAP_2")

umap1 <- cbind(umap_df, annotation_lineage=liver_milo$annotation_lineage) %>%
  ggplot(aes(UMAP_1, UMAP_2, color=as.character(annotation_lineage))) +
  geom_point(size=0.1, alpha=0.5) +
  ggrepel::geom_text_repel(data = . %>%
              group_by(annotation_lineage) %>%
              summarise(UMAP_1=mean(UMAP_1), UMAP_2=mean(UMAP_2)),
            aes(label=annotation_lineage), color="black", size=6
            ) +
  scale_color_manual(values=getPalette(colourCount)) +
  guides(color="none") +
  xlab("UMAP1") + ylab("UMAP2") +
  coord_fixed() +
  theme_classic(base_size = 22) +
  theme(axis.text = element_blank(), axis.ticks = element_blank())

umap2 <-
  cbind(umap_df, condition=as.character(liver_milo$condition)) %>%
  ggplot(aes(UMAP_1, UMAP_2, color=condition)) +
  geom_point(size=0.1, alpha=0.5) +
  scale_color_brewer(palette="Set1", name='') +
  xlab("UMAP1") + ylab("UMAP2") +
  coord_fixed() +
  guides(color='none') +
  facet_wrap(condition~., ncol=1) +
  theme_nothing(font_size = 22) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position=c(0.9,0.9),
        strip.background = element_rect(color=NA), strip.text = element_text(size=22))

nh_graph_pl <- plotNhoodGraphDA(liver_milo, milo_res, alpha = 0.1) +
  theme(legend.text = element_text(size=22), legend.title = element_text(size=24)) +
  coord_fixed()

fig4_top <- (umap1 | umap2 | nh_graph_pl) +
  plot_layout(widths = c(3,1,3))

fig4_top +
  ggsave("~/milo_output/liver_embedding.pdf", width=15, height = 10)

Next, we can check the cell types where we observe most differences between healthy and cirrhotic cells, by taking the most frequent cell type in each neighbourhood.

# Add annotation of most frequent cell type per nhood to milo results table
milo_res <- annotateNhoods(liver_milo, milo_res, "annotation_indepth")
anno_df <- data.frame(liver_milo@colData) %>%
  distinct(annotation_lineage, annotation_indepth)
milo_res <- left_join(milo_res, anno_df, by="annotation_indepth")

We first check that neighbourhoods are quite homogeneous

frac_hist <- ggplot(milo_res, aes(annotation_indepth_fraction)) +
  geom_histogram(bins=30) +
  xlab("Fraction of cells in \nmost abundant cluster") +
  ylab("# neighbourhoods") +
  theme_bw(base_size=14)

frac_hist

I can recover all the clusters where DA was detected in the original paper (see all the barplots for each lineage) and more! All in a single analysis, and without knowing where the subclusters are. Let’s bear in mind that positive logFC

Close-up on Endothelial lineage

endo_milo <- scater::runUMAP(liver_milo[,liver_milo$annotation_lineage=="Endothelia"],  dimred='PCA')
plotUMAP(endo_milo, colour_by = "annotation_indepth")

umap_df <- data.frame(reducedDim(endo_milo, "UMAP"))
colnames(umap_df) <- c("UMAP_1", "UMAP_2")

endo_umap <- cbind(umap_df, condition=endo_milo$condition) %>%
   ggplot(aes(UMAP_1, UMAP_2, color=condition)) +
  geom_point(size=0.3, alpha=0.5) +
  scale_color_brewer(palette="Set1", name='') +
  xlab("UMAP1") + ylab("UMAP2") +
  coord_fixed() +
  guides(color="none") +
  facet_wrap(condition~., ncol=1) +
  theme_nothing() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position=c(0.9,0.9),
        strip.background = element_rect(color=NA), strip.text = element_text(size=22))

liver_milo2 <- liver_milo
subset.nhoods <- str_detect(milo_res$annotation_indepth, "Endo")
reducedDim(liver_milo2, "UMAP")[colnames(endo_milo),] <- reducedDim(endo_milo, "UMAP") 

endo_gr <-
  plotNhoodGraphDA(
  liver_milo2, milo_res,
  subset.nhoods = subset.nhoods, 
  # ) =)[1:(length()-1)], 
  alpha = 0.1
  )  +
   theme(legend.text = element_text(size=22), legend.title = element_text(size=24))
  
fig4_bright1 <- ((endo_umap + endo_gr ) + 
  plot_layout(widths = c(1,2), 
                guides = "collect"
                )) 
fig4_bright1 +
  ggsave("~/milo_output/liver_endoGraph.pdf", width=9, height = 5)  

Close-up on Cholangiocytes

Filter out cells that show contamination from immune cells (expression of immune markers)

keep <- logcounts(chol_milo)["CD19",] == 0 | logcounts(chol_milo)["MS4A1",] == 0
chol_milo <- chol_milo[,keep]
chol_milo <- scater::runUMAP(chol_milo,  dimred='PCA')

plotUMAP(chol_milo, colour_by = "annotation_indepth")

umap_df <- data.frame(reducedDim(chol_milo, "UMAP"))
colnames(umap_df) <- c("UMAP_1", "UMAP_2")

chol_umap <- cbind(umap_df, condition=chol_milo$condition) %>%
   ggplot(aes(UMAP_1, UMAP_2, color=condition)) +
  geom_point(size=0.3, alpha=0.5) +
  scale_color_brewer(palette="Set1", name='') +
  xlab("UMAP1") + ylab("UMAP2") +
  coord_fixed() +
  guides(color="none") +
  facet_wrap(condition~., ncol=1) +
  theme_nothing() +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position=c(0.9,0.9),
        strip.background = element_rect(color=NA), strip.text = element_text(size=22))

chol_umap

liver_milo2 <- liver_milo
subset.nhoods <- str_detect(milo_res$annotation_indepth, "Cholangio")
reducedDim(liver_milo2, "UMAP")[colnames(chol_milo),] <- reducedDim(chol_milo, "UMAP") 

chol_gr <-
  plotNhoodGraphDA(
  liver_milo2, milo_res,
  subset.nhoods = subset.nhoods, 
  # ) =)[1:(length()-1)], 
  alpha = 0.1
  )  +
   theme(legend.text = element_text(size=22), legend.title = element_text(size=24))
  
(chol_umap + chol_gr ) + 
  plot_layout(widths = c(1,2), 
                guides = "collect"
                )

# fig4_bright1 +
#   ggsave("~/milo_output/liver_endoGraph.pdf", width=9, height = 5)  

DGE analysis

Feature selection

dec_liver <- modelGeneVar(liver_milo)
fit_liver <- metadata(dec_liver)
hvgs <- getTopHVGs(dec_liver, n=3000)

Add nhood expression to speed-up plotting of heatmaps

liver_milo <- calcNhoodExpression(liver_milo, assay = "logcounts", subset.row = hvgs)

Endothelia

unique(nhs.da.gr)

[1] 1 2 3

Visualize as volcano

marker.df %>%
  mutate(label=ifelse((adj.P.Val_2) < 0.01, GeneID, NA)) %>%
  ggplot(aes(logFC_2, -log10(adj.P.Val_2), 
             # color=highlight
             )) + 
  geom_point() +
  ggrepel::geom_text_repel(aes(label=label)) +
  xlab("logFC") + ylab("- log10(Adj. P value)") +
  theme_bw(base_size = 22)

marker.df %>%
  mutate(label=ifelse((adj.P.Val_1) < 0.01, GeneID, NA)) %>%
  ggplot(aes(logFC_1, -log10(adj.P.Val_1), 
             # color=highlight
             )) + 
  geom_point() +
  ggrepel::geom_text_repel(aes(label=label)) +
  xlab("logFC") + ylab("- log10(Adj. P value)") +
  theme_bw(base_size = 22)

NA
NA

Visualize as heatmap

(gene expression values are scaled between 0 and 1 for each gene)

GO term analysis

em_res_up

em_res_down

Cholangiocytes

x = liver_milo
da.res = milo_res
subset.row = hvgs
assay="counts"
aggregate.samples = TRUE
sample_col = "dataset"
gene.offset=TRUE

subset.nhoods = milo_res$annotation_lineage=="Cholangiocytes"
overlap=1

nhood.adj <- nhoodAdjacency(liver_milo)[as.character(unlist(nhoodIndex(liver_milo)[subset.nhoods])),as.character(unlist(nhoodIndex(liver_milo)[subset.nhoods]))]

nhood.adj[which(milo_res[subset.nhoods, "logFC"] > 0),which(milo_res[subset.nhoods, "logFC"] < 0)] <- 0
nhood.adj[which(milo_res[subset.nhoods, "logFC"] < 0),which(milo_res[subset.nhoods, "logFC"] > 0)] <- 0

nhood.adj[nhood.adj < overlap] <- 0 
g <- graph_from_adjacency_matrix(nhood.adj)
nhs.da.gr <- components(g)$membership


nhood.gr <- unique(nhs.da.gr)
# perform DGE _within_ each group of cells using the input design matrix
message(paste0("Nhoods aggregated into ", length(nhood.gr), " groups"))

Nhoods aggregated into 3 groups

fake.meta <- data.frame("CellID"=colnames(x), "Nhood.Group"=rep(NA, ncol(x)))
rownames(fake.meta) <- fake.meta$CellID

# do we want to allow cells to be members of multiple groups? This will create
# chaos for the LM as there will be a dependency structure comparing 2 different
# groups that contain overlapping cells.
# this approach means that the latter group takes precedent.
# maybe exclude the cells that fall into separate groups?
for(i in seq_along(nhood.gr)){
    nhood.x <- nhs.da.gr == nhood.gr[i]
    # get the nhoods
    nhs <- nhoods(x)
    if(!is.null(subset.nhoods)){
        nhs <- nhs[subset.nhoods]
    }

    if(!any(is.na(fake.meta[unlist(nhs[nhood.x]),]$Nhood.Group))){
        fake.meta[unlist(nhs[nhood.x]),]$Nhood.Group[!is.na(fake.meta[unlist(nhs[nhood.x]),]$Nhood.Group)] <- NA
        } else{
            fake.meta[unlist(nhs[nhood.x]),]$Nhood.Group <- nhood.gr[i]
        }
}

# only compare against the other DA neighbourhoods
x <- x[, !is.na(fake.meta$Nhood.Group)]
fake.meta <- fake.meta[!is.na(fake.meta$Nhood.Group), ]

if(!is.null(subset.row)){
    x <- x[subset.row, , drop=FALSE]
}

exprs <- assay(x, assay)

marker.list <- list()
i.contrast <- c("TestTest - TestRef") # always use contrasts for this

    # if there is only 1 group, then need to make sure that all neighbourhoods
    # are not in this group - otherwise can't do any DGE testing
    if(length(nhood.gr) == 1){
        if(sum(fake.meta$Nhood.Group == nhood.gr[1]) == nrow(fake.meta)){
            warning("All graph neighbourhoods are in the same group - cannot perform DGE testing. Returning NULL")
            return(NULL)
        }
    }
    
    ## Aggregate expression by sample
    # To avoid treating cells as independent replicates
    if (isTRUE(aggregate.samples)) {
        fake.meta[,"sample_id"] <- colData(x)[[sample_col]]
        fake.meta[,'sample_group'] <- paste(fake.meta[,"sample_id"], fake.meta[,"Nhood.Group"], sep="_")
        
        sample_gr_mat <- matrix(0, nrow=nrow(fake.meta), ncol=length(unique(fake.meta$sample_group)))
        colnames(sample_gr_mat) <- unique(fake.meta$sample_group)
        rownames(sample_gr_mat) <- rownames(fake.meta)
        
        for (s in colnames(sample_gr_mat)) {
            sample_gr_mat[which(fake.meta$sample_group == s),s] <- 1  
        }
        
        ## Summarise expression by sample
        exprs_smp <- matrix(0, nrow=nrow(exprs), ncol=ncol(sample_gr_mat))
        if (assay=='counts') {
            summFunc <- rowSums
        } else {
            summFunc <- rowMeans
        }
        
        for (i in 1:ncol(sample_gr_mat)){
            if (sum(sample_gr_mat[,i]) > 1) {
                exprs_smp[,i] <- summFunc(exprs[,which(sample_gr_mat[,i] > 0)])  
            } else {
                exprs_smp[,i] <- exprs[,which(sample_gr_mat[,i] > 0)]
            }
        }
        rownames(exprs_smp) <- rownames(exprs)
        colnames(exprs_smp) <- colnames(sample_gr_mat)
        
        smp_meta <- distinct(fake.meta, sample_group, Nhood.Group)
        rownames(smp_meta) <- smp_meta[,"sample_group"]
        
        fake.meta <- smp_meta
        exprs <- exprs_smp
    }
    
    
    for(i in seq_along(nhood.gr)){
        i.meta <- fake.meta
        i.meta$Test <- "Ref"
        i.meta$Test[fake.meta$Nhood.Group == nhood.gr[i]] <- "Test"

        if(ncol(exprs) > 1 & nrow(i.meta) > 1){
            i.design <- as.formula(" ~ 0 + Test")
            i.model <- model.matrix(i.design, data=i.meta)
            rownames(i.model) <- rownames(i.meta)
        }

        sink(file="/dev/null")
        gc()
        sink(file=NULL)

        if(assay == "logcounts"){
            i.res <- .perform_lognormal_dge(exprs, i.model, model.contrasts=i.contrast,
                                            gene.offset=gene.offset)
        } else if(assay == "counts"){
            i.res <- .perform_counts_dge(exprs, i.model, model.contrasts=i.contrast,
                                         gene.offset=gene.offset)
            colnames(i.res)[ncol(i.res)] <- "adj.P.Val"
        } else{
            warning("Assay type is not counts or logcounts - assuming (log)-normal distribution. Use these results at your peril")
            i.res <- .perform_lognormal_dge(exprs, i.model,
                                            model.contrasts=i.contrast,
                                            gene.offset=gene.offset)
        }

        i.res$adj.P.Val[is.na(i.res$adj.P.Val)] <- 1
        i.res$logFC[is.infinite(i.res$logFC)] <- 0

        i.res <- i.res[, c("logFC", "adj.P.Val")]
        colnames(i.res) <- paste(colnames(i.res), nhood.gr[i], sep="_")
        marker.list[[paste0(nhood.gr[i])]] <- i.res
    }

    marker.df.chol <- do.call(cbind.data.frame, marker.list)
    colnames(marker.df.chol) <- gsub(colnames(marker.df.chol), pattern="^[0-9]+\\.", replacement="")
    marker.df.chol$GeneID <- rownames(i.res)

table(nhs.da.gr)

nhs.da.gr
  1   2   3 
114  40   1 

Visualize as volcano

Visualize as heatmap

(gene expression values are scaled between 0 and 1 for each gene)

GO term analysis

em_res_up_chol

---

Assemble figure

fig4_bottom <- ((fig4_bleft + plot_layout()) |
      ((fig4_bright1 + plot_layout(tag_level = 'keep')) / (fig4_bbright + plot_layout())) +
      plot_layout(heights = c(1,1.6))
   ) +
  plot_layout(widths=c(1,1.2))

(fig4_top / fig4_bottom) +
  plot_layout(heights=c(1,1.8))  +
  ggsave("~/milo_output/fig4.pdf", height = 26, width = 22, useDingbats=FALSE)
  # ggsave("~/milo/ms/figures/figs/figure5.pdf", height = 26, width = 22, useDingbats=FALSE)

Assemble supplementary figure